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Abstract. A powerful method for calculating the eigenvalues of a Hamiltonian 

operator consists of converting the energy eigenvalue equation into a matrix equation by 
means of an appropriate basis set of functions. The convergence of the method can be 
greatly improved by means of a variational parameter in the basis functions determined 
by the principle of minimal sensitivity. In the case of the quartic anharmonic 
oscillator and of a symmetrical double-well potential we choose an effective oscillator 
frequency. In the case of nonsymmetrical potential we add a coordinate shift in a 
two-parameter variational calculation. The method not only gives the spectrum, but 
also an approximation to the energy eigenfunctions. Consequently it can be used to 
solve the time-dependent Schrodinger equation using the method of stationary states. 
We apply it to the time development of two different initial wave functions in the 
double- well slow roll potential. 
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1. Introduction 

We present a method for obtaining arbitrarily precise approximations to the solution of 
the time-dependent Schrodinger equation with a potential which fulfills the condition 
lim|a,|_^oo V{x) = +00 (i.e. a potential which only admits bound states). Although there 
are many examples of potentials of this kind, only a limited number of them can be 
solved exactly, the best known example being the simple harmonic oscillator (SHO), 
which is a standard topic in virtually any quantum mechanics textbook and which can 
be used to model many physical systems. 

In this paper we consider problems which cannot be solved analytically and where 
an alternative strategy must be found. Perturbation theory is the standard tool which 
is used to deal with such problems; unfortunately, the straightforward application of 
perturbation theory to some problems is not practical because the perturbation series is 
divergent. There are various methods to overcome this apparent limitation; for example 
the linear delta expansion (LDE) E] and other variants of variational perturbation 
theory (VPT) [3 IH IHl El- Loosely speaking, these techniques, although differing 
in the details, are based on the powerful idea that one can obtain a new expansion 
in some "unnatural" parameter (i.e. one not appearing in the original problem) and 
that the sequence of approximants resulting from this expansion can be made to 
converge very fast by suitably choosing a variational parameter. For example, the 
linear delta expansion works by interpolating the full Hamiltonian with the Hamiltonian 
corresponding to a soluble model, which depends on an arbitrary parameter, and by 
applying perturbation theory to it. The parameter is then determined by means of the 
principle of minimal sensitivity (PMS) |8j. Since the optimum value of the adjustable 
parameter given by the PMS depends upon the natural parameters in the Hamiltonian, 
the result corresponds to a non-perturbative result, i.e. to a non-polynomial expression 
in the natural parameters. 

Among other applications of the LDE and VPT we mention an improved Lindstedt- 
Poincare method UHl lUI; the calculation of the period of classical oscillators |T21 
[131 EI, the spectrum of a quantum potential with the WKB method (THj and the 
acceleration of the convergence of mathematical series ^H] ■ However, it was found that 
the LDE fails to give the correct long-time behavior of the wave-function in the quantum 
mechanical version of the slow-roll scenario of inflation ^7]. In successive orders it is able 
to approximate the exact time development more and more accurately, but only up to 
the time where the wave function has spread out and is beginning to contract again. The 
Hartree-Fock method does give a general qualitative picture of the time-development at 
later times, but is very far from being accurate. 

Here we propose an alternative method that consists of converting the eigenvalue 
equation into a matrix equation by taking matrix elements with respect to harmonic 
oscillator wave functions of arbitrary frequency Q, and then determining Q by some 
version of the PMS criterion. The particular criterion used here was first proposed 
and utilized in fHl- Having thus obtained the approximate energy eigenvalues and 
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eigenfunctions, one can then use the method of stationary states to calculate the time- 
dependence of the state for a given initial configuration. It turns out that this method 
is extremely accurate at even quite small orders, and has no problem with the long-time 
behaviour. 

The article is organized as follows: in Section |21 we describe the method in general 
terms and apply it to calculating the spectrum of various anharmonic oscillators; in 
Section El we use the method of stationary states to follow the time development of 
two initial wave-functions for the double-well potential that has been used in slow roll 
inflation and compare our results with those found in the literature; finally, in Section 
IHwe draw our conclusions. 



2. Energy Spectrum 

2.1. The Method 

We tackle the problem of solving the energy eigenvalue equation 

HtPn = En^n (l) 

by converting it to a matrix equation, using an orthonormal basis of wave functions of 
the quantum harmonic oscillator, depending upon an arbitrary frequency Q = a"^: 

= iV„ e-"'^'/2 ^ (2) 

where the normalization constant is = (q;/(2"' n\ y/n))^. H is the Hamiltonian for a 
particle in a one- dimensional potential that supports only bound states. 

It is necessary to truncate the infinite-dimensional matrix Hne to some finite 
dimension, say N x N, and then its eigenvalues can be calculated by simple matrix 
diagonalization. It is known that as N increases the approximation for the energy levels 
should steadily improve. This is indeed the case, but for an arbitrary Q the convergence 
may be quite slow. This leads one to seek for some criterion to choose an optimum value 
of Q. The criterion we shall adopt here, which is essentially that adopted in ^H]; is the 
principle of minimal sensitivity [S] applied to the trace of the truncated matrix. 

The rationale behind this principle is that the eigenvalues and other exact quantities 
of the problem are independent of Q but any approximate result coming from the 
diagonalization method for finite N exhibits a spurious dependence on the oscillator 
frequency. This also applies to the trace of the matrix, the sum of the eigenvalues. 
However, for finite a spurious fl dependence will emerge in the sum. A reasonable 
criterion for choosing is therefore to take it at a stationary point of 7/v = Y.n=o Hnn, 
so that this invariance is respected, at least locally. Thus we impose the PMS condition 

The reason for applying this condition to the trace is that T/v is a simple quantity to 
evaluate, and moreover it is invariant under the unitary transformation associated with 
a change of basis. Once VL is so determined, one obtains an approximation to the first 
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eigenvalues and eigenvectors of by a numerical diagonalization of the truncated NxN 
matrix. One could also contemplate applying the PMS to the determinant, which of 
course shares the property of invariance under unitary transformations, but this would 
be a much more cumbersome calculation, and could well introduce many spurious PMS 
values. 

In the following sections we will need the harmonic oscillator matrix elements of 
x^. Closed formulas have been given in Ref. which we adapt here for completeness. 

^ (2v/^)2- ^0 2^^-f'-\r - X - ky.{n - k)\{2X + ky.kl ^ ' 

for i — n = 2A, and 

^ (2v/^)2'-+i fr'o 2^--''~^+^r-X-k)\{n-k)\{2X + l + k)\k\ ^' 

for i — n = 2X + 1. These formulas assume i > n, but the matrix is symmetric. In both 
cases r must be greater than A. For all other values of p, n and i the matrix elements 
vanish. 

In addition to these we will need the matrix elements of p'^, which are given by 
1 



{p )ne = -7:" 



(6) 



2.2. The quartic anharmonic oscillator 

We take the Hamiltonian as 

i7 = ^(p2 + mV) + (?x^ (7) 

whose matrix elements in the basis of the wave functions of Eq. ^ are 

1 1 

Hnt = -{p^)ne + -m^{x'^)ne + g{x^)n£- (8) 
The trace to order is given by 

±.r„^Nin-^)^^^ii^2N'). (9) 

It turns out that 7/v has a single minimum, located at 

^PMS = 1^ - l^'^^ (10) 



-9^(1 + 2N'^) + VS (-iV^m^ + 27^2(1 + 2K 



2\2\ 2 



:ii) 



where 

-4 

The graph of this trace (divided by A^) against Q is shown in Fig. ^ exhibiting 
a global minimum. Here we have taken the parameters as m = 1, g = 1000, and 
later multiplied the eigenvalues by a factor of 2 for comparison with the results of [20] 
corresponding to (3 = 2000. By taking Q at the minimum, in accordance with PMS, we 
obtain a precise approximation to the spectrum. 
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Figure 1. Trace of the hamiltonian matrix normalized by the number of states as 
a function of the variational parameter fi. We use the parameters m = 1, g = 4000, 
corresponding to the case studied in [211 ■ Different curves represent different subspaces. 

The error in the energy of the nth. excited state, i.e. A = \E^^ — E^^^'^^l, decreases 
exponentially with N, as is shown in Fig. |21for the ground state and for the 49th excited 
state. E^^^ is the energy of the nth state obtained by considering the N x N subspace, 
whereas E^^^'^^ is the energy of the nth state corresponding to the largest subspace 
considered in this paper, which is used as reference. 

For = 100 we obtain the energy of the ground state of the Hamiltonian 
+ + 2(yfx^ as 

= 13.38844170100806193900617690280728652296098988517435666039 9 

which agrees in the first 58 (underlined) digits with the corresponding result of 
Table 1 of obtained with a different method. Much higher precision can be 
obtained by enlarging the N x N subspace: this costs little additional effort because, 
once the PMS has been applied, the Hamiltonian matrix is fully numerical and its 
eigenvalues/eigenvectors can be calculated numerically with accuracy and speed. 

It is well known that the accuracy of the eigenvalues given by the method of 
Rayleigh-Ritz decreases as the quantum number increases: this happens because the 
influence of the states which lie outside the subspace x is felt more strongly by 
the "border states" , i.e. those states which fall on the border of the selected subspace. 
In order to calculate highly excited states without enlarging the subspace too much 
one can center the subspace around that particular state and apply the procedure 
mentioned above. This has been done in Fig. El where we have plotted the ratio 
An = {E^^ — E^ ^^^) / E^ considering a square subspace of elements Hij centred 
around the n = 100 state. E^^^^ is obtained using the analytical formula for the 
energy of the quantum anharmonic oscillator obtained in PH]- The flattening of the 
error around n for the largest subspaces considered (81 x 81 and 161 x 161) signals that 
the numerical results obtained with the present method have reached the precision of 
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Figure 2. Error in tlie energy of the ground and 49*'* excited states as a function of 
n and for subspaces of different dimensions. We use the parameters m = I, g = 4000, 
corresponding to the case studied in |2(Jj . 
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Figure 3. A„ = [E^l^^ — E^'n^^^^) / E^^^^^ as a function of the dimension of the 
subspace. 



the WKB formula. We conclude that the present method can be used to obtain the 
energies and wave functions of arbitrarily high excited states with the desired accuracy. 
As expected, the error A„ is maximal for the "border states" , 

When the centred subspace is restricted to one element = 1 the method reduces 
to a well-known simple variational calculation. 

2.3. The double-well potential 

The same method can be applied with very little change to the Hamiltonian which has 
been used in the consideration of slow roll inflation in the early universe, namely 

iJ = ip2 ^ \{x^ - a^f/2A + const. 
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= -p^ - -m^x^ + (12) 

with = and g = A/24. 

The parameters will be taken as a = 5 and A = 0.01, as these have been used in 
previous work on the subject (see, for example, [HI |221 EB]) 

All that is needed in this case is to reverse the sign of in the formulas given in 
the previous subsection. We obtain similar accuracy for the eigenvalues with very little 
effort. These will be used in Section 3 to give the time development of a given initial 
configuration. 



2.4- General anharmonic potentials 

We now consider general anharmonic potentials of the form: 

V{x) = Y.^^^^' (13) 
i=o 

where the coefficients Kj define a polynomial of order M (even). We require that nj^ > 
to ensure that only bound states are permitted. 
The shifted potential 

V{x + a) = f: {x + ay = E E (l) x^-'a' , (14) 

and the original one have the same spectrum. Therefore we can choose a to be a 
variational parameter and apply the PMS as in the case of the adjustable oscillator 
frequency. 

As an example, let us consider the potential: 

V{x) = 11- 118s - 44x2 + 80x^ + 16x^ . (15) 

Since V{x) is strongly asymmetric we expect that decomposing it with respect to a basis 
of SHO wave functions centred at the origin will not be the best choice. We therefore 
translate the potential[2I] by an arbitrary quantity —a and impose the PMS on Q and 
a simultaneously. That is, we impose the two PMS conditions 

Using a 10 X 10 subspace we obtain the optimal values cr = —3.889 and Q = 31.179. 
Correspondingly we find the energy of the ground state to be: 

= - 1229.11605104 5 (17) 

where the underlined digits are correct. Using a 40 x 40 subspace we find a = —3.583 
and Q = 27.431, obtaining 

= - 1229.11605104600459705899 2 (18) 

where the underlined digits are correct. Similar results hold for the excited states. 
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We notice that apMs does not correspond to the value for which V{x) has a 
minimum, i.e. x ~ —3.979. In fact (JpMS tends to increase as the dimension of the 
subspace is increased, as a resuh of the influence of the highly excited states, which are 
less localized. 

It is clear that such a scheme can be applied with limited effort to a general 
anharmonic potential of arbitrary order: in fact, a suitable choice of a allows one to 
reach high precision with a limited number of terms. Clearly, since a general potential 
can be always expanded in a Taylor series around a point, this implies that our method 
can be easily applied to calculate portions of the spectrum of a potential with arbitrary 
precision. 

In Fig. El we have plotted the first 200 states of the spectrum of Eq. (fT3)) . using our 
method with subspaces of 100 x 100 and 200 x 200 and 500 x 500, together with the 
first-order WKB estimate for comparison. 



> 




Figure 4. Potential V{x) = 11 - 118a; - Ux^ + SOa;^ 





Figure 5. Spectrum of the potential of Eq. H15() . as a function of the quantum number 
n. The pluses, triangles and crosses correspond to the eigenvalues obtained using our 
method with subspaces of 100 x 100, 200 x 200 and 500 x 500 respectively. The first- 
order WKB estimate is also shown for the higher states. 
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3. Time development 

3.1. The Method of Stationary States 

If solutions of the energy eigenvalue equation (P) are known for a given Hamiltonian H 
then the time-dependent Schrodinger equation {h = 1) 

i—^ = H'^ . (19) 
dt 

can be solved by the method of stationary states. Namely, if the initial wave-function 
at t = can be expanded as 



^(x,t = 0) = 5]a„V'„(a;), (20) 

?i=0 

its value at a later time is given by 

oo 

m{x,t) = Y.ane-'^-' i:n{x) (21) 

?i=0 

The method of Section |21 gives us an approximation not only to the spectrum, but 
also to the energy eigenf unctions, namely 

N-l 

i^nix) = dnk (Pkix) . (22) 

A;=0 

where dnk denotes the kth component of the nth eigenvector of the truncated 
Hamiltonian matrix. 

Similarly, the initial wave function can also be expressed as a truncated expansion 
in terms of the 4>n{x)'- 

N-1 

^nix,t = 0)=Y.CnMx). (23) 

fc=0 

By comparison with Eq. ()22|) we see that 

N-1 

an = I] Q (d"^)fa . (24) 

£=0 

These coefficients can now be used in Eq. ()2H) to obtain an approximation to the 
wave function \E'(x, t) at any later time t. 

3.2. Slow roll inflation 

Here we use the Hamiltonian of Eq. (jl2|) with two different initial configurations. 
3.2.1. Centred Gaussian 

The initial wave function, used in previous studies of slow-roll infiation, is given by 

nx.t = ^) = [^)"'\—''\ 
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Our task is to find the coefficients c„ in Eq. (j22I)- 
By orthonormality, 

/ 771 \ 1/4 r+00 

(l)lix)^ix,t = 0)dx = — / e"^ ^ (25) 



we obtain 



where 



,27r 

where = m/2 + a^. By means of the change of variable y = ax 

Cn = — 7^ / e Hn{y)dy . (26) 
and the expansion of Hn{y), namely 

= g(-l)'^;^^2"-,»- . (27) 

(/— \ n— 2fc+l 

In fact c„ = unless n is even n = 2i, and then 

The coefficients a„ are now given by Eq. 

For comparison with previous work, we use our time- dependent wave function to 
calculate (x^), given by 

= j m*{x,t) x^ ^!{x,t) dx . (30) 

In terms of the ipn{x) this is 

=Y^ ala, e-*""^* 1 V'n(a;) ^,(x) rfx (31) 

where = En — E^. Using Eq. this becomes 

{^^) = E < a,e-^--* E dnkd,j{x\, . (32) 

The result for (a;^)^/^ is plotted in Fig. IHl As can be seen, the result is vastly superior 
to Hartree-Fock, and on this scale cannot be distinguished from that obtained using 
Fourier transform methods. The ratio of different orders is shown in Fig. [7| 

3.2.2. Shifted Gaussian 

In this case we consider an initial wave function 

m{x,t = 0) = 
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Figure 6. [x^)^/'^ versus t for the slow roll potential of Eq. (|12() . Various orders 
of our method compared with the Hartree-Fock method used in 23 and the Fourier 
transform method of |24| . 




Figure 7. Ratio of the results of our method for the slow roll potential: orders 10 and 
20 versus 40. 



representing a particle localized around x = Xo at t = 0. As before, we obtain the 
coefficients c„ by orthonormality: 



= Q)dx 

/ U \ /■+00 /32 



(x- 



2/3^ 



where now = /i/2 + a^. With a change of variable to y = [3{x — 



(3 \2tx) 



(33) 



(34) 
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We expand the Hermite polynomial, to obtain 

'"/^l n' / ax V^'' 



fe=0 



(n-2A;)!fc! V-^y ^ 2/3 



Therefore 



/? V27ry 

-£g'-^) \,,(.-2.-.), (i) (t) 



where 



ir,- = y^e^y'/^dy = 2^^+^^'^ \l + (-1)^1 T{j + 1/2) . (36) 
Finally we obtain 



n-2k 



" /5 V2vr; to ,t^o ^ ^ A;!(2j)!(n-2fc-2j)! 

(\ n~2k / \ n—2k—2j 

I) (t) (^)'"'^0--V2). (37, 

Figure |H1 shows the expectation value (x) as a function of time for four different 
values of /i. As expected, as /i increases, the frequency of oscillation between the two 
wells increases as well. The plots in Fig. |H1 correspond to the values g = 1/2400, 
m = l/(2v^) and Xo = 5. 




Figure 8. {x) with four different clioices of /i and witli Xq = 5 for the slow roll 
potential with a shifted Gaussian. 
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4. Conclusions 

We have shown that the matrix method, combined with the principle of minimal 
sensitivity, is a powerful tool for finding the spectrum of arbitrary polynomial potentials 
having only bound states. The original method proposed in ^H] has been generalized in 
two ways: for the calculation of higher levels, the matrix can be centred around a higher 
levels of the SHO, and for noneven potentials the introduction of a shift parameter 
improves the accuracy of the method for a given dimension A^. This aspect of the 
method could be used for a general, non-polynomial potential, in conjunction with its 
Taylor expansion about a given point. 

As a by-product of the calculation of the spectrum, the method also provides 
approximations to the energy eigenfunctions. Knowing both the energies and their 
eigenfunctions we are then able to implement the method of stationary states to track 
the time development of a given initial wave function with good accuracy and for long 
time scales, as we have demonstrated for the slow-roll potential. 

P.A. acknowledges support of Conacyt grant no. C01-40633/A-1 and of Alvarez- 
Buylla fund of the University of Colima. A. A. acknowledges support from Conacyt 
grant no. 44950 and PROMEP. 
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